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Axially symmetric relativistic MHD simulations 
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On the origin of torus and jet-like features 
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Abstract. The structure and the evolution of Pulsar Wind Nebulae (PWNe) are studied by means of two-dimensional axisym- 
metric relativistic magnetohydrodynamic (RMHD) simulations. After the first imaging of the Crab Nebula with Chandra, a 
growing number of objects has been found to show in the X-rays spatial features such as rings and jets, that clearly cannot 
be accounted for within the standard framework of one-dimensional semi-analytical models. The most promising explanation 
suggested so far is based on the combined effects of the latitude dependence of the pulsar wind energy flux, shaping the wind 
termination shock and naturally providing a higher equatorial emission, and of the wind magnetization, likely responsible for 
the jet collimation by hoop stresses downstream of the shock. This scenario is investigated here by following the evolution of a 
PWN interacting with the confining Supernova Remnant (SNR), from the free expansion to the beginning of the reverberation 
phase. Our results confirm the oblate shape of the wind termination shock and the formation of a polar jet with supersonic 
velocities (v ~ 0.5 - 0.7c) for high enough values of the equatorial wind magnetization parameter (cr > 0.01). 
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1. Introduction 

Pulsar Wind Nebulae (PWNe, or plerions) arise from the con- 
finement of pulsar winds by the surrounding medium, usually 
an expanding Supernova Remnant (SNR). The relativistic mag- 
netized pulsar wind is slowed down to non-relativistic veloci- 
ties at a termination shock, where the magnetic field is com- 
pressed and the bulk energy of the outflow is converted into 
heat and acceleration of particles. These then give rise to the 
synchrotron and Inverse Compton emission observed from ple- 
rions in a very wide range of frequencies, extending from radio 
wavelengths to X-rays and even y-rays. 

The best studied plerion is the Crab Nebula, whose emis- 
sion has been extensively investigated in all frequency bands 
and for which most models have been proposed. New light on 
the spatial structure of the Crab Nebula emission at high fre- 
quencies has been shed by observations made with the Chandra 
X-ray satellite (Weisskopf et al. 2000 1, which, thanks to the 
unprecedented spatial resolution, has revealed a number of in- 
triguing features in the inner part of the nebula (see also Hester 
et al. 1995 2002 1. The new details highlighted strengthen the 
view of the Crab Nebula as an axisymmetric object. In what 



is thought to be the equatorial plane of the pulsar rotation, 
Chandra observations show the presence of a bright ring of 
emission, lying at a much closer distance to the pulsar than 



the already identified X-ray torus (e.g. Hester et al. 1995). The 
most puzzling discovery, however, is probably the presence of 
two opposite jet-like features oriented along an axis perpen- 
dicular to the plane of the torus and emerging from the very 
close vicinity of the pulsar. Similar features have been observed 
also in a number of other objects, namely around the Vela pul- 
sar (Helfand et al. 1207771 Pavlov et al. 1277011 . PSR 1509-58 
(Gaensler et al. 2002) and in the supernova remnants GO. 9+01 
(Gaensler et al. l200Tl and G54. 1+0.3 (Lu et al. l2002l 

While the presence of a X-ray bright torus may be at least 
qualitatively explained within the framework of standard 1- 
D RMHD models (Kennel & Coroniti 173531 KC84 hereafter; 
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Emmering & Chevalier 1987|), if we further assume that ei- 
ther the energy flux emerging from the pulsar or the termi- 
nation shock dissipation efficiency is higher at low latitudes 
around the equator, the presence of jets that seem to emanate 
directly from the pulsar poses severe theoretical problems in 
its interpretation (Lyubarsky & Eichler 2001 ), given the diffi- 
culties at explaining self-collimation of ultra-relativistic flows. 
A recent suggestion for an answer to this puzzle (Bogovalov 
& Khangoulian 2002a 2002b Lyubarsky 2002 ) is that the jets 
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are actually originating downstream of the pulsar wind termi- 
nation shock, where the flow is only mildly or non-relativistic. 
If this is the case, the fact that they are observed starting from 
a much closer distance from the pulsar than where the shock 
in the equatorial plane is thought to be, has to be interpreted 
assuming that the given degree of anisotropy in the energy flow 
from the pulsar also causes the shock front to be highly non- 
spherical in shape, much closer to the pulsar along the rotation 
axis than in the equatorial plane. Moreover, even if the pulsar 
wind is weakly magnetized just upstream of the termination 
shock, the magnetic field inside the plerion can become as high 
as to reach equipartition. Therefore, collimation of the down- 
stream flow may be easily achieved there by magnetic hoop 
stresses (Lyubarsky 2002, Khangoulian & Bogovalov I2003> . 
resulting in plasma compression toward the axis and eventu- 
ally in a polar jet-like outflow. 

Thanks to the recent progress in numerical relativistic fluid 
dynamics and MHD (see Del Zanna & Bucciantini 2002 Del 
Zanna et al. 2003 and references therein), we are now able to 
start a more quantitative investigation of this problem by means 
of computer simulations. Our aim is to clarify whether a given 
latitude dependence of the pulsar wind energy flux may actu- 
ally explain the jet-torus morphology observed at X-ray fre- 
quencies for the Crab Nebula and other plerions, and, if this 
is the case, what are the conclusions that one may infer on 
the structure and magnetization of the unshocked pulsar wind. 
Here we present the results of a first series of long-term 2-D ax- 
isymmetric RMHD simulations, from which some general con- 
clusions on the physical mechanisms at work and useful scal- 
ings may already be derived (see Amato et al. 2003 for prelim- 
inary results). A similar numerical investigation has been re- 
cently carried out (Komissarov & Lyubarsky 2003 KL03 here- 
after), confirming the basic physical picture as viable for ex- 
plaining the main observational features, as strongly suggested 
also by the close resemblance, at least at a qualitative level, be- 
tween the map of simulated emission and Chandra images of 
the Crab Nebula. 

The paper structure is as follows. In Sect.[2]the pulsar wind 
model adopted at large distances from the light cylinder is 
sketched. In Sect. [3] the numerical details of the simulations 
and the initial conditions are reported. Sect. |4] deals with the 
results of the simulations, split in three sub-sections for con- 
venience. Finally the results are summarized in Sect. [5] where 
conclusions are drawn for this preliminary work. 

2. Pulsar wind model and pre-shock conditions 

The key point of all attempts at interpreting the torus and the 
jet-like features observed in plerions as arising post-shock is 
that the energy flux in the unshocked pulsar wind should de- 
pend on latitude as sin 2 8 (Bogovalov & Khangoulian 2002a 
Lyubarsky 2002). This angular dependence is related to the 
structure of the residual purely toroidal magnetic field, which, 
far enough from the pulsar, is expected to depend on the po- 
lar angle as sin 6* (split monopole models: e.g. Michel [19731 
Bogovalov [1999). We further assume that along the way to 
the termination shock the residual Poynting energy flux is al- 
most entirely converted into particle energy flux, as in clas- 



sic models, preserving the overall angular dependence. We do 
not address here the fundamental problem of how the acceler- 
ation of the outflow and the conversion of magnetic to parti- 
cle energy may occur, the so-called cr paradox: see however 
Vlahakis (20041 and references therein for self-similar 2-D 
relativistic MHD models, though many other mechanisms in- 
volving waves, reconnection or kinetic plasma processes have 
been proposed (see e.g. Arons 2003 for a review and references 
therein). Therefore, in the present work we will consider as an 
initial condition an axisymmetric cold (p <k pc 2 ) relativisti- 
cally expanding pulsar wind (v r ~ c), with a small ratio be- 
tween electromagnetic and particle energy fluxes just upstream 
of the termination shock, thus extending the standard KC84 
picture to 2-D. 

Following the previous analytical studies cited above, let us 
then introduce the latitude dependence of the energy flux from 
the pulsar as a dependence on 6 of the wind Lorentz factor y, 
namely: 

y(#)=y [a + (l-ff)sin 2 #], (1) 

where the subscript indicates quantities in the equatorial 
plane, and a < 1 is a parameter controlling the ratio be- 
tween the Lorentz factor at the pole and that at the equator. We 
then assume the streamlines to be radial upstream of the shock 
and the mass flux to be isotropic (Bogovalov & Khangoulian 
I2()()2a> : 

with ro the (arbitrary) distance from the pulsar on the equatorial 
plane at which the value po is assigned. As it was shown in the 
above cited work, these assumptions, in the case of a weakly 
magnetized outflow, naturally give rise to an oblate shape of 
the termination shock front, and imply the existence, around 
the axis, of a colder, denser sector of the downstream plasma. 

As anticipated, the toroidal magnetic field B = B$ is de- 
fined as 

B(r,6) = B — sinfl. (3) 
r 

It is important to notice that the adopted pulsar wind model 
leads to a ratio between Poynting and kinetic energy fluxes that 
does depend on latitude but not on distance from the pulsar. 
Therefore, it is convenient to take its equatorial (maximum) 
value cr as the independent parameter controlling the wind 
magnetization (thus the same quantity used by KC84) and let 
Z?o derive from it as Bo = (^PQC 2 y^cr) l l 2 . The total wind en- 
ergy flux F ^ c(pc 2 y 2 + B 2 /An) may then be written as 

F(r, 9) = po c 3 7o (jf [a + d - a + cr) sin 2 6], (4) 

and the wind anisotropy in energy flux is, at any distance, given 
by the value F(r,Q)/F(r,7r/2) = a/(l + cr). Finally, the value 
of po is determined from the total pulsar spin-down luminosity, 
supposed to be constant in time for simplicity (see Bucciantini 
et al. 2003 , 2004 for long-term spherically symmetric simula- 
tions in the case of a decaying luminosity in time). 
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A similar wind model and pre-shock conditions were 
adopted by KL03, though the assumed mass flux was not 
isotropic but had the same latitude dependence as the energy 
flux (y — jo — 10 and F(8) ~ sin 8 2 ), and the toroidal magnetic 
field a different shape: while in our case we take a field with 
its maximum at the equator, allowing direct comparison with 
standard 1 -D models and the usual definition of cr, in the cited 
work the field reaches the maximum at intermediate latitudes 
and then goes smoothly to zero at the equator (see Sect. l4.3l for 
a detailed discussion and comparison). In spite of these differ- 
ences, the overall latitude dependence of the wind energy flux, 
which is the quantity that really matters in shaping the termi- 
nation shock, is in both cases the one predicted in the cited 
analytical studies. 

3. Simulation setup 

3.1. Numerical settings 

The problem of the interaction of the pulsar wind sketched 
above with the surrounding medium is here addressed numer- 
ically by performing 2-D axisymmetric simulations. The nu- 
merical tool employed is the shock-capturing code for rela- 
tivistic MHD developed by Del Zanna et al. ( HUOll EOTM 
The scheme is particularly simple and efficient, since complex 
Riemann solvers based on characteristic waves are avoided in 
favor of central-type component-wise techniques: the solver is 
defined by the two fastest local magnetosonic speeds and spa- 
tial reconstruction at cell boundaries is achieved by using ENO- 
type interpolating polynomials. However, due to the extremely 
small time-steps required in the present simulations and be- 
cause of the high Lorentz factors involved (numerical oscilla- 
tions are most dangerous in the ultra-relativistic regime), third 
order reconstruction is avoided and simpler second order lim- 
ited reconstruction is employed. Time integration is achieved 
by means of a two-stage Runge-Kutta TVD algorithm, using a 
CFL number of 0.5. 

The spatial grid is defined by spherical coordinates with 
400 cells in the radial direction and 100 cells in the polar an- 
gle 8. The physical domain of the simulation ranges in radius 
from r m ; n = 0.05 to r max = 20 light-years (we assume a unit 
length r = 1 light-year, from now on all lengths will be ex- 
pressed in light-years, time intervals in years and velocities in 
units of c, if not explicitly stated otherwise), and a logarith- 
mic stretching (dr ~ r) is imposed to better resolve the inner 
region. Note that, in order to resolve within a few computa- 
tional cells the contact discontinuity between the lighter rel- 
ativistic plasma and the heavy SNR ejecta (typical jumps of 
order 10 , see below), artificial compression is required, espe- 
cially at large radii where resolution is necessarily lower. This 
may in principle amplify spurious noise above the threshold of 
dissipation by numerical viscosity, thus producing unphysical 
results or even code crashing. We have verified that the cho- 
sen grid and scheme settings are a good compromise between 
resolution and stability. 

Stationary input for all quantities is imposed at r m j n , where 
the super-Alfvenic wind is blowing from, while zeroth order 
extrapolation is assumed at the outer boundary. The domain in 



8 is the first quadrant (0, n/2), with reflecting conditions for vg 
and B at the polar axis to enforce axisymmetry, and on vg alone 
on the equatorial plane. Note that all quantities are defined in 
our code at cell centers, which never coincide with the symme- 
try axis 8-0 where singularities may occur. By doing so, the 
ghost cells technique is well suited for all non-cartesian prob- 
lems and we never find carbuncle-like effects (see the jet propa- 
gation tests in cylindrical geometry in Del Zanna & Bucciantini 
12001 Del Zanna et al. l2003l . 

Notice that under these particular settings, only 5 (out of 
8) RMHD variables need to be evolved in time. Moreover, the 
magnetic field is forced to be in the azimuthal direction alone, 
and thus always perpendicular to the velocity vector, by the 
assumed geometrical symmetries. In this case, all the specific 
methods to treat the divergence-free constraint are obviously of 
no use, and the magnetic field may be evolved as an ordinary 
fluid variable. Other important simplifications occur in the al- 
gorithm for calculating the local magnetosonic speeds and in 
that for deriving the primitive variables from the set of con- 
servative variables (see Del Zanna et al. 2003 1. In order to 
speed up calculations, the pre-shock quantities are not updated 
in time and the global time-step is defined at the inner termina- 
tion shock radius (on the axis), which is moving outward thus 
having the effect of accelerating the simulation. 

3.2. Initial conditions and simulation parameters 

The initial conditions are given as follows. The pulsar wind, 
modeled as in Sect. |2] is set up within an arbitrary radius of 
0.2. Once the two free parameters a and cr are provided, the 
wind is uniquely determined by the equatorial Lorentz factor yo 
and total wind luminosity Lq. Appropriate values of the above 
parameters come by fits of the high-energy emission (for the 
Crab Nebula) based on spherically symmetric MHD (KC84) 
and kinetic (Gallant & Arons 1994 1 models of the wind-nebula 
interaction, namely yo > 10 6 and Lq - 5 x 10 38 erg s" 1 . 

From a numerical point of view, however, such a high 
Lorentz factor is well beyond the capabilities of existing rela- 
tivistic MHD codes. Here we assume yo = 100 (notice that ours 
are the multidimensional RMHD simulations with the highest 
Lorentz factor available in the literature so far), and an aver- 
aged wind luminosity of 5xl0 39 erg s , in order to speedup the 
evolution. The resulting rest mass densities will be of course 
unrealistically high, since we basically have poyjj = const - We 
deem that this should not a problem even in multidimensional 
flows, provided y » 1 in the wind region at all latitudes (as 
in our case). This property has been nicely demonstrated in 1- 
D by KC84, where there is no explicit dependence on y and 
p taken separately, while all quantities depend on the wind lu- 
minosity and on the magnetization parameter cr. In the present 
2-D case, even if the termination shock is not circular in shape 
(thus the wind is not always normal to it and KC84 analysis do 
not apply at all latitudes), we still find that the exact value of 
yo does not make any difference in the nebular flow, provided 
all other parameters (except obviously po) are kept the same. 
In particular we have performed runs with lower values of the 
equatorial Lorentz factor, down to yo = 20 (here with a small 
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wind anisotropy, a = 0.5, to preserve the condition y » 1 at 
the poles too), and we did not find appreciable differences in 
the results. 

In all simulations we assume an anisotropy parameter of 
a = 0.1, thus the energy flux along the polar axis will be ten 
times less than at the equator, while density will be ten times 
greater. We decided not to choose smaller values of a since 
we wanted to be sure that during the evolution the termination 
shock (TS hereafter) would eventually move away from the in- 
ner boundary at all latitudes, otherwise outflow conditions are 
not appropriate and may lead to unphysical situations, typically 
near the polar axis. The magnetization parameter cr is instead 
varied in the different simulations, from 0.003 up to 0.1. 

Around the pulsar wind region, the (spherical) expansion 
of the cold, unmagnetized supernova (SN) ejecta is simulated 
by setting up a region of high constant density with a radially 
increasing velocity profile v = v e jr/r e j, appropriate for a self- 
similar expansion. Here we take r e i = r^ — 1, while p e j and v e j 
are respectively determined by imposing 
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where we take M ej = 3M G = 6x 10 g, and 



Jo 2 



p e j v 2 4nr 2 dr, 



(5) 



(6) 



with E e j = 10 erg. The velocity at the outer boundary of the 
ejecta is then v e j a 7500 km s~ , corresponding to an age of 
the SNR of r e jlv e j « 40 years, while the velocity at the contact 
discontinuity between the ejecta and the relativistic material 
(CD hereafter) is v e j w 1500 km s . Farther out, between r e j 
and r max , ISM conditions are imposed, that is a uniform, static, 
unmagnetized background with pism = 10~ 24 g cm -3 Pism ~ 
10~ 12 dyne cm . For similar settings see van der Swaluw et al. 
(2001 1, Blondin et al. SMTH . Bucciantini et al. 120051120031 . 

Finally, for simplicity we adopt here a uniform value of 4/3 
for the adiabatic index, appropriate for a relativistic plasma. 
Radially symmetric simulations of the PWN - SNR interaction 
with a variable adiabatic index (5/3 in the ejecta and ISM re- 
gions) were performed by Bucciantini et al. ( 2003 1, to which 
we refer for a discussion of the complications implied. 

Note that in KL03 the ISM is absent and ejecta expand- 
ing with a velocity of 5000 km s -1 are set everywhere beyond 
the wind region. This allows to speed up the evolution of the 
TS, though important processes due to the interaction with the 
external ISM (namely the reverberation phase, see below) are 
completely neglected. 

4. Simulation results and discussion 

4.1. Overall PWN structure and evolution: comparison 
with analytical models 

Before studying the formation of the peculiar jet-torus structure 
seen in Chandra images, let us investigate the overall plerion 
properties and its evolution in time, comparing the results with 
analytical models, whenever possible. 
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Fig. 1. The time evolution of the PWN boundaries, the TS and 
CD radii, for 8 = n/2 and 6 = (symbols as indicated on the 
plot), in the cr = 0.003 case. The 2-D results are compared 
with a spherically symmetric case (solid curves) set up with 
the parameters corresponding to the equatorial outflow of the 
2-D case. The self-similar 1-D hydro solution (dashed curves) 
is also shown for comparison. 

The PWN evolution is followed up to t = 1000 for four 
cases with different magnetization: cr = 0.003, cr = 0.01, 
cr = 0.03, and cr = 0.1. After a short (a few years) transient 
stage during which, after the nebula is first formed, the reverse 
shock propagates backward, both the wind termination shock 
and the contact discontinuity (the latter separates the nebula 
from the swept up shell of ejecta) move outward. In Fig.^the 
evolution of the PWN boundaries for 8 = tt/2 and 8 = is 
plotted against time, in the cr = 0.003 case. For a comparison, 
also the corresponding 1-D spherically symmetric evolution is 
shown, together with the fits expected for (hydrodynamical) 
self-similar models of PWN interacting with freely expanding 
SN ejecta (see Bucciantini et al. 2004 and references therein). 
At later times (f » 500 in this case) the expected self-similar 
expansion is slowed down because of the interaction with the 
reverse shock produced by the motion of the SNR in the sur- 
rounding ISM. This is the beginning of the so called reverber- 
ation phase (see Bucciantini et al. 20031, here occurring rather 
early because of the high spin-down luminosity adopted. 

As expected, the PWN inner boundary (the termination 
shock, TS hereafter) is farther from the pulsar at the equator 
than at the pole, while the opposite occurs at the outer bound- 
ary (the contact discontinuity, CD hereafter). The former effect 
is due to the assumed wind energy flux anisotropy which pro- 
duces the oblate shape of the TS. The latter effect is due, in- 
stead, to the pinching by the PWN magnetic field (Begelman 
& Li 119921 van der Swaluw 2003 1. Similar results are found 
also for more magnetized winds, although for higher values 
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Evolution of TS shape 
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Fig. 2. The time evolution of the TS radius in the case cr = 
0.003, in the cylindrical coordinates R and Z (solid line). 
Together with the shock shape resulting from the simulation, 
the expression in Eq. [8] is plotted for comparison as a dashed 
curve. Before reverberation starts (f =s 500), the ratio between 
the polar and equatorial radii is also well reproduced. 
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Fig. 3. The flow structure around the TS. The background 2-D 
gray-scale plot refers to the velocity magnitude. The arrows in- 
dicate the streamlines. Labels refer to: A) ultrarelativistic wind 
region; B) subsonic equatorial outflow; C) equatorial super- 
sonic funnel; D) super-fastmagnetosonic shocked outflow; a) 
termination shock front; b) rim shock; c) fastmagnetosonic sur- 
face. 



of cr the TS evolution is slower and its expansion can begin 
only when 2-D effects (vortexes) remove the cr < vqdIc con- 
straint for quasi-stationary radial MHD flows (KC84; see also 
Bucciantini et al. 2004 for discussions about this constraint). 



Evolution of PWN elongation 
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Fig. 4. The time evolution of the PWN elongation, that is the 
ratio of the CD radii at the pole and at the equator, for the var- 
ious values of the wind magnetization parameter. The elonga- 
tion starts decreasing in time when the CD reaches the SNR- 
ISM reverse shock along the polar axis. 

Let us discuss in greater detail the above results. The time 
evolution of the TS shape is shown in Fig.[2]for the case cr = 
0.003. If the downstream total pressure were constant, the TS 
profile at a given time would be simply defined by the condition 



F(r, ff) cos 6 = const, 



(7) 



where 6 is the angle between the shock normal and the ra- 
dial direction. We have verified that the approximation of con- 
stant downstream pressure is reasonable only within an angle e 
around the equator, with e ~ 20° for a — 0.1 as in the present 
simulations. At these latitudes one also has 5*0, leading to 
the following approximate expression for the TS profile: 



Rts(0) - Rts(tt/2) 



a + ( l - a + cr) sin (6) 



l+o- 



1/2 



(8) 



The latter expression is plotted against the numerical solution 
for the shock front in Fig.|5]for each time. The deviations from 
the predicted shape at intermediate latitudes are mostly due to 
pressure variations in the post-shock region: the pressure con- 
stancy is verified only within an order of magnitude due to the 
magnetic pinching effect and the presence of supersonic flows. 
Due to this fact, no improvement is observed if one plots the ex- 
act solution of Eq.Qinstead of the expression in Eq.|S] Another 
thing to notice from the profiles shown in Fig.|5]is that the evo- 
lution of the shock shape appears to be self-similar. We have 
checked that this result holds also for all the other values of cr 
considered (as long as the PWN evolution remains in the free 
expansion stage). 

The detailed structure of the flow in the vicinities of the 
TS is shown in Fig. [3 where its complexity is apparent. Here 
we show the different regimes of the post-shock flow: in region 
D, due to the obliquity of the TS, the speed remains super- 
fastmagnetosonic, until the plasma crosses the rim shock (la- 
beled in figure as b: see also KL03), and it is finally slowed 
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down to sub-fastmagnetosonic, yet still supersonic, speeds in 
the funnel C. The TS front between the equator and the rim 
shock latitude e (see above) is almost perpendicular to the wind 
(here is where Eq. is a good approximation), so that the lat- 
ter is directly slowed down to subsonic speeds (region B). 

In Fig.Hthe PWN elongation, defined as R C d(0)/R cd (7t/2), 
is shown as a function of time for all the different values of 
the wind magnetization, thus extending the analysis by van der 
Swaluw (2003 1 to the relativistic case and to higher magne- 
tizations. We confirm the result that the elongation increases 
with time and with cr. The growth with time is limited to the 
free expansion phase, before the interaction with the reverse 
shock from the ISM. When this interaction starts, it has two 
consequences: the reduction in time of the elongation, eventu- 
ally toward unity, and the saturation of the elongation increase 
with cr. We reiterate that this interaction begins early in time 
in our simulations because of the high spin-down luminosity 
assumed. 

The situation for the low cr case is displayed in Fig. [5] 
where 2-D plots of various quantities are shown for t = 400 
and t = 1000 (the magnetization is defined here as the ratio 
between the magnetic and thermal pressure). 

At early times the PWN is clearly elongated along the sym- 
metry axis, since from the density and pressure plots we can see 
that the interaction with the outer shell (SNR blast wave and 
reverse shock) has not started yet, except for the pole where it 
is just beginning. Notice that not only the elongation, but also 
the expected shape for a PWN expanding in fast moving SN 
ejecta is, at least qualitatively, reproduced (see Fig. |2 case c, 
in the paper by Begelman & Li 1992 1. The other prediction 
in the cited work, namely the dependence of the total pressure 
on the cylindrical radius alone inside the PWN, is reproduced 
quite well, as the corresponding contours are nearly parallel 
to the vertical axis. As mentioned previously, it is only in the 
close vicinities of the TS that major deviations, including dis- 
continuities, occur. Notice, however, that such overall pressure 
distribution is supposed to hold only approximately, since the 
assumed wind energy flux is strongly anisotropic and vortexes 
of non-negligible speed may form. The presence of these vor- 
texes is also the cause of a dragging of high density material 
from the ejecta, that starts close to the pole and later extends 
to the entire nebula (see the first plot in Fig. [3}. However, this 
fact does not seem to affect the overall dynamics, at least for 
low cr values. From the third plot we can have an estimate of 
the relative importance of magnetic and thermal effects, whose 
ratio increases with the cylindrical radius, again as expected. 

At much later times, namely t = 1000, the PWN has ex- 
panded so much that the interaction with the reverse shock in 
the ejecta occurs at all latitudes. The dependence of the ther- 
mal and magnetic pressures on the cylindrical radius still holds 
approximately, but the shape of the PWN is no more elon- 
gated, since in spite of the pinching effect, which is still at 
work, the external boundaries are now defined by the spherical 
SNR shell. Due to the high value of the sound speed (actually 
fast magnetosonic speed) inside the PWN, waves and distur- 
bances created at the interaction interface between the PWN 
and the SNR are rapidly transmitted and distributed through- 



out the whole PWN, back to the TS that changes in shape and 
slows down its evolution (see Fig.[0. 

In the above discussion reference was made to the simu- 
lation with a — 0.1 and cr = 0.003 (the latter is the standard 
value deduced for the Crab Nebula from 1-D radial models, 
see KC84). Let us now briefly discuss how the above results 
depend on these parameters. The energy flux anisotropy is cru- 
cial for determining the TS shape, see Fig.[2]and Eq. (jSJl, though 
runs with different values of a show that the overall morphol- 
ogy and evolution of the PWN are basically unchanged. More 
important is the value of the wind magnetization parameter cr, 
which determines the PWN elongation and, as a side effect, 
also the time at which the interaction with the reverse shock 
begins (see Fig. @}. In the following section we will discuss 
in greater detail the role played by the magnetization in the 
physics of the PWN and in particular in the formation of the 
polar jets. 

4.2. Formation of jet-like features 

Let us now investigate how the flow pattern inside the PWN 
is affected by the nebula magnetization. In Fig. [6] we show the 
speed magnitude and the streamlines for increasing values of cr, 
namely 0.003, 0.01 and 0.03, at the same time t = 400. The jet 
is basically absent in the low magnetization case, where only a 
subsonic flow (v < 0.1) is observed along the polar axis. In the 
intermediate case the polar outflow starts to be more collimated 
and its speed increases to supersonic velocities. Finally, in the 
high magnetization case, a strong, well collimated jet is clearly 
apparent. It has supersonic speed reaching values as high as 
v* 0.7 -0.8. 

The presence of the polar jet seems to be directly corre- 
lated to the flow pattern in the rest of the nebula. In the first 
two cases an equatorial flow with speeds v ~ 0.5 is present, 
together with large scale vortexes at higher latitudes. The fast 
equatorial flow is entirely due to 2-D effects, since the stream- 
lines bend toward the equator after crossing the toroidal TS 
(Bogovalov & Khangoulian 2002a I. The vortexes occur when 
this equatorial flow hits the expanding CD boundary and a cir- 
culating back-flow is created. This pattern is present in hydro 
simulations as well with the unmagnetized case being very sim- 
ilar to the cr = 0.003 case in Fig. [6] With increasing value of 
the wind magnetization the equatorial outflow, for this choice 
of Poynting flux dependence on latitude, is progressively sup- 
pressed. For cr = 0.03 the flow in the equatorial plane is limited 
to the close vicinities of the TS. A situation very similar to the 
latter is found in the highest magnetization case we considered, 
cr = 0.1, which is not displayed in the figure. 

As mentioned above, a polar subsonic outflow is present 
even in the hydro simulations. The reason for this is that large 
scale vortexes eventually reach the high latitude regions near 
the axis, compress the plasma there and drive a polar flow. In 
addition to this process, as the cr of the wind increases, the 
magnetic hoop stresses cause the formation of the jet-like fea- 
ture with super-fast-magnetosonic speeds: the tension of the 
toroidal magnetic field, amplified downstream of the TS, di- 
verts the equatorial flow toward the polar axis, well before the 
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large scale vortexes due to interaction with the CD outer bound- 
ary are even formed. 

This mechanism has been theoretically predicted by 
Lyubarsky (2002 1, analysed by Khangoulian & Bogovalov 
(2003 1 and finally confirmed numerically by KL03. In particu- 
lar, in the last cited work the authors show that significant polar 
velocities are found when the effective wind magnetization pa- 
rameter is cr ss 0.01 or higher. In spite of the different settings 
(especially the shape of the pulsar wind magnetic field: see next 
section) it is interesting to notice that we basically confirm here 
a similar value for the threshold between cases in which the po- 
lar outflow is well developed and supersonic and cases in which 
it is not. 

Let us consider the case cr = 0.03 and follow more closely 
what happens in the inner regions around the TS. In Fig. Q 
we plot the velocity field, magnetization and density at two 
different times: t = 200 (upper row) and t = 300 (lower 
row). By following the streamlines in the first panel we can 
understand the launching mechanism. Due to the highly non- 
spherical shape of the termination shock, narrow channels of 
supersonic post-shock flow form along its boundary, and then 



converge toward the equatorial plane, where they merge in the 
high speed (v * 0.7) equatorial beams present also in low-cr 
and hydro simulations. However, this flow is no longer able to 
reach the outer boundary of the PWN. After a few TS radii, the 
hoop stresses due to the high magnetization in the inner part of 
the nebula (even beyond equipartition: see the second plot) are 
so strong as to inhibit this flow completely, thus a back-flow 
and vortexes form there. When this back-flow reaches the polar 
axis, part of it directly escapes along the axis, another part of 
it circulates back into the equatorial flow, and yet another part 
heats up the plasma in the cusp region just on top of the TS, 
driving again a polar outflow. The final result is then the com- 
plete suppression of the equatorial flow at a small radius and 
the formation of a high velocity (v ~ 0.7) polar jet. 

At t — 300, additional effects have occurred. The TS has 
moved farther out, while the equatorial outflow region has 
started shrinking, with the result of enhancing the jet driving 
mechanism even more. The reason for this unexpected behavior 
can be appreciated by the density maps. Fingers of cold, dense 
material protruding from the ejecta through the CD crunch the 
highly magnetized region of the PWN. Rather than the result 
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Fig. 6. Flow magnitude (gray scale images) and streamlines at time t = 400 for three values of the wind magnetization parameter 
cr. The jet starts to form for cr = 0.01 and it is very well developed for higher values. 



of some instability (Rayleigh-Taylor and/or Kelvin-Helmoltz) 
this density enhancement coming from the external equatorial 
regions is rather due to the large scale circulation induced by 
the jet interaction with the slowly expanding ejecta. The polar 
flow is suddenly diverted along the CD all the way down to the 
equatorial region, where it goes back toward the center eventu- 
ally dragging with it the dense material from the ejecta (faint 
vortexes with v = 0.1 — 0.2 are present also in the last plot of 
Fig.EJ. 

While the inner nebula region is very well resolved, due 
to the logarithmic radial stretching of the computational grid, 
one may wonder whether the outer large scale vortexes and 
the dragging of the dense material from the CD are indeed 
numerically converged features, since resolution becomes cor- 
respondingly poor at large distances. In order to check this 
point we have performed a long-term run with double resolu- 
tion (800 x 200). The results show that the overall structure is 
maintained, especially in the inner region (corresponding to the 
X-ray most luminous region). However, the higher resolution 
allows the development of smaller scale vortexes and instabili- 
ties, so that the dragging is more efficient and saturation occurs 
earlier, although the global behaviour is preserved. 

To summarize, the pulsar wind anisotropy is responsible for 
the non-spherical shape of the TS, which makes the material 
flow along its boundary to form eventually the high velocity 
outflow in the equatorial plane. If the wind magnetization pa- 
rameter cr is high enough, equipartition is reached in the PWN 
equatorial region very near the TS and this flow can be com- 
pletely suppressed by hoop stresses and diverted toward the po- 
lar axis, where part of it will drive the super-fastmagnetosonic 
jet. When this polar outflow starts interacting with the slowly 
expanding ejecta, the material escapes along the CD down to 
the equator, where high density cool plasma is dragged by the 
large scale vortexes toward the center causing the crushing of 
the inner part of the nebula. 



4.3. Dependence on the wind magnetic field shape 

When defining the wind properties we have assumed a sin 9 lat- 
itude dependence for the toroidal field, thus the magnetization 
is taken to be highest at the equator. However, since the toroidal 
magnetic field is the residual of the global field at large dis- 
tances from the pulsar, around the equatorial plane it is bound 
to change sign. Even assuming the simplest case, that is the 
split monopole model (Michel Tl973> for a perfectly aligned 
magnetic rotator, at large enough distances from the light cylin- 
der r » c/Q {Q. is the pulsar angular velocity) the toroidal field 
eventually becomes 



and since B, changes sign at 9 — n/2 also B$ is bound to do so. 
If the pulsar dipolar field is not exactly aligned with the rota- 
tion axis, the so-called striped wind on the equator will present 
a sequence of fieldlines of alternating sign. If magnetic dis- 
sipation (reconnection) is efficient, the effective toroidal field 
at large distances will be approximately zero below a certain 
small angle (depending on the obliquity of the rotator) around 
the equator (e.g. Coroniti fl990> . 

In our analysis we have considered an approximately 
aligned rotator, so that the striped wind region is so small to 
be dynamically unimportant and the assumption B ~ sin 9 re- 
mains valid even for 9 — > n/2 (the change of sign is superflu- 
ous, since B appears only squared in the dynamical equations). 
As mentioned in Sect.|2] in KL03 the opposite scenario is con- 
sidered: the striped wind region is very large, such as to extend 
to intermediate latitudes affecting the overall field shape. 

To understand what differences are to be expected, we have 
performed also a series of simulations with a magnetic field 
given as 

B(r, 9) = B — sinfl tanh[b(n/2 - &)], (10) 
r 

where Bq is the same parameter defined in Sect. |2 in terms of 
cr, and b is a free numerical parameter defining the width of the 
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Fig. 7. The mechanism responsible for jet formation. Flow, magnetization and density maps are shown for the cr = 0.03 case at 
two different times: t = 200 (upper row) and t = 300 (lower row). 



striped wind region. In Fig.[S]we show the angular dependence 
of B for various values of b: b — > oo corresponds to the case 
discussed in all the previous sections, with B oc sin 9, b = 10 
is the case that will be considered in this section, while the 
field used in KL03 basically coincides with that in Eq. dlOt i for 
b * 0.7. 



In Fig. [5] we show the same quantities as in Fig. at time 
t = 300, again for the case cr = 0.03 but using the magnetic 
field of Eq. fllOi with b — 10. The presence of an unmagnetized 
region around the equator (where the magnetic field was the 
highest in the previous cases) even as narrow as in the present 
case, is sufficient to change the picture completely. Now the TS 
is able to move further out at the equator. Around the narrow 
channel to which the equatorial outflow is confined, the mag- 
netization grows rapidly and the hoop stresses are effective. As 
a consequence, collimation of a polar outflow still occurs, al- 
though the fraction of plasma now involved is less than when 
the magnetic field was non-vanishing at the equator. In particu- 
Fig. 8. ^The latitude dependence of the toroidal field, as defined ^ (see Fig J3] for compar i son ) the high speed funnel, instead of 
by Eq.GS) (only the field above the equator is displayed). In being completely diverted t0 the axiSi splits into an equatorial 
this section the value b - 10 is employed. outflow and in a backflow toward the axis. The corresponding 
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streamlines now lie very close to the TS, and correspondingly 
also the polar jet forms much closer to the origin. 

5. Conclusions 

In this paper the structure and evolution of PWNe interacting 
with the surrounding SN ejecta has been analyzed by means of 
relativistic MHD axisymmetric simulations. Our main goal has 
been here the investigation of the mechanism originating the 
polar jets which are observed in X-rays in a growing number 
of PWNe. 

The most recent and promising analytical studies 
(Bogovalov & Khangoulian 2002b Lyubarsky 2002 1 start with 
the assumption that the pulsar wind is highly anisotropic, with 
a much larger energy flux at the equator than at the pole. Such 
a wind, interacting with the expanding SN ejecta produces a 
hot magnetized bubble with a torus-jet structure, as observed. 
The polar jet is originated by the magnetic hoop stresses in the 
PWN that, for high enough values of the wind magnetization 
parameter cr, divert part of the equatorial flow toward the axis, 
collimating and accelerating it. The first numerical RMHD 
simulations (Amato et al. 2003 KL03) confirm this scenario, 
and the simulated synchrotron emission map in the latter cited 
work strikingly resembles, at least qualitatively, the X-ray 
images of the Crab Nebula. 

Here we have made an effort to improve on those prelimi- 
nary simulations. The equatorial relativistic wind has a Lorentz 
factor as high as y — 100, the magnetization parameter cr 
goes from 0.003 up to 0.1, which leads to magnetic fields in 
the PWN well beyond equipartition. The magnetic field shape 
assumed in the wind is that proper for an aligned or weakly 
oblique rotator, while in KL03 the assumed field is far from 
the sin 6 dependence, with a very broad region of low magneti- 
zation around the equator. The evolution is followed up to the 
beginning of the reverberation phase and comparison with pre- 
vious models is made whenever possible. 

The results show that the predicted self-similar evolution of 
the external PWN boundary is well reproduced at low magne- 



tizations, and before reverberation effects begin. The expected 
elongation of the PWN due to magnetic pinching is also recov- 
ered, and we show how this effect increases in time and with cr. 
The elongation of the nebula appears to be independent on the 
wind anisotropy, measured in our model by the parameter a. 

In the inner part of the PWN the termination shock as- 
sumes an oblate shape, as expected for an anisotropic wind 
energy flux. An equatorial supersonic flow is produced in a 
complex shock structure. At intermediate latitudes, where the 
TS is highly oblique, the downstream flow is still supersonic. 
The plasma moves along the front gradually focusing toward 
the equator. This pattern holds in hydrodynamical simulations 
as well, since it is due only to the wind flux anisotropy. The 
equatorial flow that is driven by this focusing mechanism is su- 
personic, with typical velocities of v « 0.5 - 0.7c, consistent 
with the values inferred for motion in the equatorial plane of 
the Crab Nebula (Hester et al. l2002> . 

When the wind magnetization is high enough (cr > 0.01 in 
our simulations) the magnetic hoop stresses in the equatorial 
plane are so strong as to suppress completely this flow after 
a few termination shock radii. The plasma is then completely 
diverted toward the axis, where the magnetic compression fi- 
nally drives a supersonic polar outflow with velocities that are 
once again in agreement with the values observed in the Crab 
and Vela PWNe (v * 0.3 - 0.7c, see Hester et al. 120021 Pavlov 
et al. 2003 1. At later times the interaction of the polar jets with 
the contact discontinuity density jump may cause additional ef- 
fects. The flow circulates all the way along the CD from high 
latitudes down to the equatorial plane. Here, these large scale 
vortexes drag some dense material from the ejecta toward the 
origin, with the effect of confining the equatorial outflow to the 
very inner parts of the nebula. Notice that the circulation is the 
opposite (from the equator to the pole) in the hydro and low cr 
cases. The value of cr that distinguishes the two regimes that we 
call highly and lowly magnetizatized depends on the expansion 
velocity of the CD: for nebulae expanding at a larger rate we 
expect the transition to occur for higher wind magnetizations. 
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Finally, for the high cr cases, the flow pattern strongly de- 
pends on the Poynting flux distribution in the wind. In particu- 
lar, if in a narrow region around the equator the magnetic field 
vanishes, an additional vortex appears in the circulation pattern. 
The equatorial supersonic flow is never suppressed completely: 
it reaches the CD and then it circulates back toward the axis. 
The polar jet is still present and drives a vortex circulating at 
higher latitudes in the opposite direction than the former. 

In conclusion, axisymmetric relativistic MHD simulations 
of the interaction of pulsar winds with expanding SNRs are 
able to reproduce at least qualitatively most of the structures 
seen in X-ray images: the overall toroidal structure of the ple- 
rion, the supersonic motions in the equatorial plane and, if the 
wind magnetization is high enough, also the presence of po- 
lar jets with supersonic velocities. A more detailed study, with 
a larger sampling of the parameter space would be necessary 
for a quantitative comparison with the observations. An inter- 
esting perspective that the present study suggests is the infer- 
ence of the wind magnetic field structure from direct compar- 
ison between simulated synchrotron maps and X-ray observa- 
tions. Preliminary results are encouraging, although we prefer 
to leave quantitative comparisons as future work. 

Some of the observed features are anyway impossible to 
reproduce within the present axisymmetric RMHD framework. 
Emission structures like knots and sprites might well be related 
to non-ideal effects, like magnetic reconnection and dissipa- 
tion, that are non-trivial to deal with. Remaining within the 
RMHD approximation, a full 3-D setting would allow to deal 
with some non-axisymmetric instabilities that may play an im- 
portant role in PWNe: let us mention as an example the kink 
instability, which is likely to be at the origin of the bending of 
the jet (observed in both the Crab and Vela nebulae). We also 
leave the study of these important physical processes as future 
work. 
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